%% This is a temporary copy of fungi_fluocell_detect.m to experiment with colourpsace basis changing
%% input the information of the image, i.e. the path, name, and format
function eigenvectordilation(file_name)
%     p = 'C:/Users/wang-lab/Desktop/data/pathology/0622/fungi/';
%     p = 'C:\Users\wang-lab\Desktop\data\pathology\0109_2015\negative_slide_tile\';
    p = 'C:\Users\wang-lab\Desktop\data\pathology\0622\fungi\';
    fungi_intensity_th = 100;
    % Consider increasing the minimal area. 
    min_area = 500; % minimal area for fungi in image

    im = imread([p,file_name]);
    % im = uint8(im);
    [~,name,ext] = fileparts(file_name);
    figure;imagesc(im);
    
    %% reverse the image before background subtraction
    rev_im = 255 * uint8(ones(size(im))) - im;
    
    %% subtract background on counter staining
    file_bg = strcat(p,'output/background_', name,'.mat');
    [bw, poly] = get_background(im, file_bg);
    figure; imshow(rev_im);
    im_revonly = rev_im;
    hold on; plot(poly(:,1), poly(:,2),'b','LineWidth',1.5);
    for i = 1:3,
        bg_mean = compute_average_value(rev_im(:,:,i),bw);
        % Background subtraction followed by converting RGB to Index
        rev_im(:,:,i) = uint8(rev_im(:,:,i) - bg_mean);
    end;
    im_prep = sum(rev_im, 3);
    
    %% find the global threshold of the image and the fungi area
    level = my_graythresh(im_prep);
    fg_bw = im2bw(im_prep, level) & (im_prep > fungi_intensity_th);
    fg_bw = bwareaopen(fg_bw, min_area);% define a minimal area for fungi and get rid of staining dots
    fg_bwlabel = bwlabel(fg_bw);
    num_items = max(max(fg_bwlabel));
    fg_bw = my_WatershedSeg(fg_bw);
    fg_bd = bwboundaries(fg_bw, 'noholes');
    stats = regionprops(fg_bw, 'Orientation');
    orientationv = -[stats.Orientation]';
    cropped_masks = cell(1, num_items);
    cropped_images = cell(1, num_items);
    cropped_images_max_min = cell(1, num_items);
    %% get cropping masks
    for i = 1:length(num_items)
        [cropped_masks{i}, cropped_images{i}] = cropper(rev_im, fg_bwlabel, fg_bd{i});
        temp = cropped_images{i};
        temp = cropped_masks{i}.*temp(:,:,3); % we just want blue channel
        min = 255; 
        dim = size(temp);
        max_val = max(max(temp));
        max_find = true; % FIGURE OUT A METHOD WITHOUT FOR LOOPS
        max_coords = 0;
        min_coords = 0;
        for x = 1:dim(1)
            for y = 1:dim(2)
                if(max_find && temp(x, y) == max_val)
                    max_coords = [x y];
                    max_find = false;
                elseif(temp(x, y) ~= 0 && temp(x, y) < min) 
                    min = temp(x, y);
                    min_coords = [x y];
                end;
            end;
        end;
        cropped_images_max_min{i} = [max_coords min_coords];
    end;
    %% eigenvalue dilation
    lambda = 5;
    dilation_matrices = cell(1, length(num_items));
    for i = 1:length(num_items)
        cur = cropped_images{i};
        coords = cropped_images_max_min{i};
        x = [cur(coords(1), coords(2), 1); cur(coords(1), coords(2), 2); cur(coords(1), coords(2), 3)]; 
        x_prime = [cur(coords(3), coords(4), 1); cur(coords(3), coords(4), 2); cur(coords(3), coords(4), 3)];
        x_sub = x - x_prime;
        v1 = double(x_sub);
        v1 = v1 + [1; 0; 1];
        v2 = double(x_sub);
        v2 = v2 + [0; 1; 0];
        P = orth(double([(x-x_prime) v1 v2]));
        D = diag([lambda 1 1]);
        dilation_matrices{i} = P*D*(P\eye(3));
    end;
return;